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ABSTRACT 

A  computer  program  is  presented  for  the  small  strain  analysis  of 
plane  structures  in  the  strain  hardening  elastic-plastic  range.  The 
finite  element  displacement  method  is  used  to  perform  the  linear 
analyses  in  the  iterative  scheme.  Bar  and  constant  strain  isotropic 
plane  stress  triangles  are  available  for  use  in  constructing  idealiza¬ 
tions.  The  use  of  ten  different  sets  of  material  properties,  three 
different  material  laws,  and  incremental  proportional  loading  are 
available  as  options.  Good  correlation  is  shown  with  available  test 
data  and  theoretical  solutions. 

The  distribution  of  this  abstract  is  unlimited. 
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SECTION  I 
INTRODUCTION 

This  report  describes  a  computer  program  for  the  stress  analysis 
of  plane  structures  in  the  elastic -plastic  range  by  the  finite  element 
method.  The  program  can  handle  bar  and  triangular  plate  elements  so 
that  it  is  applicable  to  trusses  and  to  the  analysis  of  in-plane  stresses 
in  reinforced  plates.  The  material  behavior  is  assumed  to  be  isotropic 
and  the  user  has  a  choice  of  three  types  of  stress -strain  laws  and  ten 
different  materials. 

A  numerical  step  by  step  procedure  for  obtaining  solutions  which 
satisfy  the  requirements  of  the  incremental  theory  of  plasticity  for 
materials  which  obey  the  Mises  yield  condition  and  the  associated  flow 
rule  is  used  in  the  program.  At  each  step  in  the  solution,  an  iterative 
procedure  is  used  to  find  the  correct  values  of  the  strain  increments. 
Changes  in  plastic  strain  are  accounted  for  by  the  addition  of  fictitious 
plastic  forces  to  the  actual  loading  on  the  structure  in  such  a  way  that 
the  deflections  of  the  structure  under  the  modified  loading  with 
assumed  elastic  behavior  are  equal  to  the  actual  deflections.  A 
modified  form  of  the  computer  program  given  in  Reference  1  is  used 
to  obtain  elastic  solutions. 

The  finite  element  method  originally  developed  for  elastic  stress 
analysis  was  extended  to  apply  to  inelastic  material  behavior  by  Padlog 
(Reference  2)  et  al.  Further  studies  of  the  use  of  the  method  in 
elastic -plastic  problems  have  been  made,  for  example,  at  MIT 
(Reference  3)  and  the  California  Institute  of  Technology  (Reference  4). 
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The  method  has  been  extended  to  apply  to  anisotropic  materials  by 
Jensen  (Reference  5).  Marcal  and  King  (Reference  7)  have  applied 
the  method  to  problems  of  plane  stress  and  strain  and  to  axisymmetric 
problems.  However,  the  computer  programs  necessary  for  the  ap¬ 
plication  of  the  method  have  not  been  published.  The  purpose  of  this 
report  is  to  make  such  a  program  available. 

i  r ..  ■  Vi  . .  vi '  . 

The  following  section  of  the  report  presents  a  detailed  explanation 
of  the  method  of  analysis.  The  application  of  the  program  is  demon¬ 
strated  by  examples  in  Section  HI.  Directions  for  the  use  of  the 
computer  program  are  given  in  Section  IV.  Listing  of  the  program 
and  sample  data  sets  are  given  in  the  Appendix. 

A , 

al  invert 
•  +  >  ’ 
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SECTION  n 

METHOD  OF  ANALYSIS 

The  analysis  procedure  used  in  the  computer  program  is  described 
here.  The  method,  first  used  by  Padlog  for  the  solution  of  problems 
involving  plastic  flow  and  creep,  is  given  here  in  a  slightly  modified 
form.  The  well  known  formulas  for  the  stiffness  of  bars  and  triangular 
plate  elements  are  first  presented.  Then  the  step  by  step  iterative 
procedure  used  for  the  solution  of  problems  in  which  plastic  flow 
occurs  is  described. 

1.  ELEMENT  PROPERTIES 

The  two  types  of  elements  to  be  considered  in  this  analysis  (the  bar 
and  the  triangular  plate)  are  shown  in  Figure  1.  The  coordinates  of 
the  end  points  of  the  bar  and  the  vertices  of  the  triangle  are  referred 
to  a  fixed  coordinate  system  in  the  plane.  The  cartesian  components 
of  the  nodal  displacements  for  each  of  these  elements  comprise  the 
element  displacement  vector  X.  The  ordering  of  the  displacement 
components  is  shown  in  Figure  1.  The  total  element  strains  designated 
by  the  vector  c  can  be  expressed  in  terms  of  the  nodal  displacements 
by  an  equation  of  the  form 

e  =  BX  (  I ) 

The  stresses  are  related  to  the  elastic  strains  c  e  =  «  —  by 

Hooke's  law 

O’  =  C€®  (2) 
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x'  =  I  ui  »v  i »  U2i  v2 


Bar  Element 


Figure  1.  Bar  and  Plate  Elements 
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The  nodal  forces  F  corresponding  to  given  displacements  X  are  found 
by  the  principle  of  virtual  work.  That  is 

X'F  =  f  tVdV  (3) 

V 

where  X  and  e  are  the  virtual  displacement  vector  and  the  corres¬ 
ponding  strain  vector,  respectively,  and  the  integration  is  carried  out 
over  the  volume  of  the  element.  Prime  superscripts  denote  transposed 
matrices.  Using  Equations  1  and  2  in  Equation  3  gives  the  results 

X  ’  [  F  -  /  B'CBdVX  +  /b'C  d  V«P  |  =0  (4) 

V  v  J 

Or,  since  the  elements  of  X  are  arbitrary,  then 

F  +  Fp  =  kX  (5) 

where 

k  =  J  B*CBdV  (6) 

V 

is  the  element  stiffness  matrix  and 

Fp  =  D  €  p  (7) 

is  the  vector  of  plastic  forces  in  which 

D  =  j  B*  C  d  V  (8  ) 

V 

The  definitions  of  the  nodal  force  and  displacement  vectors  and  of  the 
matrices  C,  B,  D,  and  k  for  bars  and  triangular  plate  elements  are: 

Bar  Element 

x‘  =  [  u,  ,  v,  ,  u2  ,  v2] 

F '  =  [  F  ,  F  ,  F  ,  F  1 
l  xi  *  yi  X2  *  y2  J 

8  =  l2  (  ~ *21  ’  ~  y2l  *  *  21  ’  *21  ] 


(9) 

(10) 
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in  which  L  is  the  length  of  the  bar  and 


X_  r  X  _  —  X  ,  V_  =  V  —  y  etc 
21  2  1*  '21  f2  '| 

C  =  E  —  Youngs  modulus 


TSmb  aojccjrr  r..h  0  =  [  -  X  ,  -  y  ,  x  ,y| 

L  l  21  *  *21  ’  21  ’  *21 J 

Wtfere  A  is  the  cross-sectional  area  of  the  bar 

r 

T»~>  - :  ~P 


’  _*  AE 

m  C  .  =--Vs  in  a  Mjiurc k  =  ■  - 

ir  ~zii  ^3 


in  which 

I'i.LI  rnn- 


k 

-  k . 


'Hke  i Diva  .  '  ^ii 

’Mux  '  ^  i-ur..'.: 

Tfciangiilar  Plate  Element 


21 


x  y 
21  21 


“k. 

k 


*21  y2l 


'21 


^  *V.  *U2*V2*U3*V3] 


:&tfhefe 


h  *21  V3I 


(12) 

(13) 

(14) 


(15) 


(16) 


(17) 


ffesd'g:  nx 

F‘  =  l 

:rc?r  _ 

v?  zt  r  •:  • 

F  ,  F  ,  F 

xi  *  y  i  ’  x 2  * 

Fy2  ’ 

F  ,  F  1 

X3  ’  y3 . 

1 

(18) 

ismi 

S •  % ■  r,,  1 

(19) 

57  ^  1c  fiGjj.'na*  | 

CT,  •  ■  Tx,  ] 

(20) 

-»J2  0 

»3. 

O 

-*2. 

- 

o 

B  =  — 

aSrratr*--  °Si  t| 

£.  3 

VP  °  *32 

0 

~  *  31 

0 

*21 

(21) 

M-ir  • 

~  *32  ~ V32 

”*31 

V3. 

*  21 

"V2. 

*3.  *2. 


(22) 
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The  absolute 


in  which  v  - 


D 


where  = 


where 


value  of  h  equals  twice  the  area  of  the  triangular  element. 

I 

v 

\-ur 


C  = 


V 

I 


0 

0 


0 

0  (1-1/ 

)/2 

Poisson's  ratio 

V32 

-  i/y 

7  32 

c( 

X 

32 

*/X32 

X  32 

-c. 

y32 

E  t  Ihl 

^31 

^y3. 

-c, 

X  31 

2h  (1  -  i/2  ) 

“  X  31 

C, 

y  3. 

- 

y2. 

~  21 

c. 

X  2 1 

vx 

21 

X  21 

-c 

1 

.  y  2. 

(l-i/)/2  and  t  is  the  element  thickness. 

E  t  - 


k  = 


2lhl  (I  -v  ) 


(23) 


(  24) 


(25) 


k.. 


(26) 
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k  1 1 

= 

y32  •  *  12  "  _l/X32y32  “  ^1 

X32  y  32 

k  13 

= 

~y3iy32_  ^lX3lX32’  *14 

=  l/X3l  y  3  2  +  Cl  y  3 1  X  3  2 

T,5 

= 

y2l  y32  +  Cl  *21  *32  ’  *  16 

'  "X2.  *32  '  Cl  y  2  1  *  32 

k  22 

= 

X  32  +  Ciy32  1  *23  "  Vy3l  X  32 

+  CIX  31  y  3  2 

*24 

= 

_X3I  X  32  ~^iy3l  y  3  2  ’  *23 

~  J/y2|  X  32  ~  ^1  X  2 1  y  32 

r26 

= 

X  21  X32  +  e.y2.  y  3  2  •  "*33  = 

y32.  +C.  X3.  * 

*34 

= 

VXZ\  y3 |  y 3 1  X 3 |  ’  *35 

=  _y2.  X  3  1  ~  C 1  X  2 1  X  31 

T36 

= 

VXZX  y  3  1  +  C.y2.  X  3  1  ’  T44 

*  c,v|,  ■ 

F45 

= 

I/y2.  X  31  +  C  1  X  21  y  3 1  •  F46 

_  *21  *31  ^  1  y 21  y3l 

*55 

= 

y  2 1  +  ^lX2l  ’  **36  “  ~VXZ\ 

y 2 1  ”  C  1  X  2 1  y2. 

*66 

X||  +C.  y2. 

2.  ELASTIC -PLASTIC  ANALYSIS 

In  the  elastic  range  of  material  behavior  the  equilibrium  equations 
for  a  structure  composed  of  plate  and  bar  elements  of  the  type  considered 
here  can  be  written 

F  =  KX  (28) 

where  the  force  and  displacement  vectors  now  have  as  their  components 
the  cartesian  components  of  force  and  displacement  at  all  the  nodes 
and  K  is  the  assembled  stiffness  matrix  for  the  whole  structure.  The 
solution  of  Equation  28  for  the  unknown  displacement  is  given 
symbolically  by 

X  =  K"1  F  (29) 
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The  displacements  known,  the  element  strains  can  be  obtained  from 
Equation  1  and  the  stresses  from  Equation  2.  However,  when  the 
stresses  reach  the  intensity  required  to  cause  plastic  flow,  it  becomes 
necessary  to  determine  the  increments  of  plastic  strain  caused  by  the 
load  increment.  The  material  is  assumed  to  obey  the  Mises  yield 
condition  and  the  associated  flow  rule.  For  plane  stress  the  following 
equations  apply 


(  cr2  -  cr  cr  +  cr2  +  3  t2  1 
\  X  x  y  y  xy  ) 


t> 

-'l 

II 

7f(HP 

A*p  = 

X 

2  <T  x 

A,P  = 

A?-  (20-,  • 

y 

2cr  y 

A/xPy  = 

,  afp 

3  ^  rxy 

1,2  =  H(«P> 


2  .1/2 


(30) 


(31) 


(32) 


where  cr  and  e  P  are  the  effective  stress  and  the  effective  plastic 
strain,  respectively,  and  where  HU  /  is  the  stress-plastic  strain 
relation  for  uniaxial  stress. 


If  it  is  assumed  that  the  response  of  the  structure  to  the  removal 
of  a  load  increment  will  be  completely  elastic  then  Equation  28  can  be 
modified  to  account  for  plastic  flow  as  follows 

K  X  =  F  +  Fp  (33) 

where  X  and  F  are  the  displacement  and  load  after  the  application  of 
the  increment  and  FP  is  the  vector  of  plastic  forces  corresponding  to 
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the  plastic  strains.  The  plastic  strain  increments  caused  by  the  incre¬ 
ment  of  load  must  satisfy  Equations  32  and  for  an  element  undergoing 
plastic  flow  the  stresses  must  satisfy  the  yield  condition  (Equation  30). 

The  following  step  by  step  iterative  method  is  used  to  obtain 
solutions: 

1.  An  increment  is  given  to  the  applied  loads. 

2.  New  values  of  displacement  are  found  from  Equation  33  using 
the  current  values  of  the  plastic  forces  (these  will  be  zero  for 
the  first  step). 

3.  The  displacements  are  used  to  compute  total  strains,  elastic 
strains,  stresses,  and  the  effective  stress. 

4.  If  the  new  value  of  the  effective  stress  is  greater  than  the 
largest  previous  value,  the  element  is  plastic  and  the 
effective  stress  is  used  to  determine  a  new  value  of  the 
effective  strain. 

5.  Plastic  strain  increments  computed  from  Equations  32  are 
added  to  the  current  values  of  the  plastic  strain  and  new 
values  of  the  plastic  forces  are  calculated. 

6.  If  the  increment  in  effective  plastic  strain  is  sufficiently  small 
the  iteration  is  complete  and  a  return  to  step  1  is  made,  if  not 
a  return  is  made  to  step  2  and  a  new  cycle  begun. 
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This  procedure  is  applied  to  each  of  the  elements  and  the  decision  to 
start  a  new  step  (apply  a  load  increment)  is  based  on  the  largest  plastic 
strain  increment  found  among  all  the  elements. 

An  important  feature  of  the  method  is  the  way  in  which  the 
effective  plastic  strain  is  computed  from  the  new  value  of  the  ef¬ 
fective  stress  at  each  iteration.  If  the  inverse  of  Equation  30  is  used 

_  D  _ 

to  give  €  as  a  function  of  a  the  solution  may  become  unstable.  This 
becomes  obvious  when  one  considers  the  case  of  the  elastic,  perfectly 
plastic  material  for  which  the  inverse  of  the  function  H(7P)  does  not 
exist.  To  avoid  this  difficulty  the  "constant  strain"  method  of 
Reference  2  is  used.  In  this  method  the  total  strain  e*  is  taken  equal 
to  the  sum  of  the  value  of  ?  p  computed  in  the  previous  iteration 
and  a/  E  . 


The  stress -strain  law  can  be  written  in  the  form 


f- 


or 


H(i-p) 


+  7P 


(34) 


The  new  value  of  e  p  can  be  found  from  Equation  34  without  difficulty. 


The  criterion  used  in  step  6  of  the  iterative  procedure  given  above, 
to  decide  whether  the  plastic  strains  have  been  determined  with  suf¬ 
ficient  accuracy,  is  the  size  of  the  ratio  of  the  increment  in  effective 
plastic  strain  to  F/E.  This  ratio  is  a  measure  of  the  difference 
between  the  ordinates  to  the  theoretical  stress  strain  curve  and  the 
curve  that  is  actually  being  used  at  that  step  in  the  calculation. 
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3.  STRESS-STRAIN  LAWS 


The  following  three  types  of  stress-strain  laws  are  available  for 
use  in  the  computer  program.  Each  of  them  is  a  three  parameter  law. 


Type  1  -  Ramberg-Osgood  Law 


in  which 


—  +  3(Ti  /  <r  \n 
E  7E  U,  ) 


E  —  Young's  modulus 

oj  —  secant  yield  stress  (stress  at  which  the 
secant  modulus  =  0.  7E) 

n  —  shape  factor 


Type  2  -  Goldberg-Richard  Law 


in  which 

E 

°i 

n 


<T 


E  €  , 


1  -L 

E«f 

n  “ 

1  “ 

_ 

—  Young's  modulus 

—  maximum  stress 

—  shape  factor 


Type  3  -  Bilinear  Law 


cr 


<T 


Ee 


t 


for  a  <  a 

y 

for  cr  >  cry 
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in  which 

E  —  Young's  modulus 

ay  —  yield  stress 

E^  —  Slope  of  the  plastic  portion  of  the 
stress-strain  curve 

To  reduce  computing  time  a  linearized  form  of  the  Ramberg-Osgood 
law  is  used  in  the  program.  This  law  is  fitted  by  a  series  of  straight 
line  segments  which  match  the  actual  curve  at  100  points  in  the  interval 
0  <  <  20  cr,  /  E  .  If  values  of  outside  this  range  are  en¬ 

countered  the  exact  formula  is  used. 
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SECTION  HI 
EXAMPLES 


1.  NONLINEAR  TRUSS 

To  illustrate  the  use  of  the  program  a  solution  for  the  member 
stresses  and  the  tip  deflection  of  the  cantilever  truss  shown  in 
Figure  2  is  obtained.  A  stress-strain  relation  of  the  Ramberg-Osgood 
type  is  assumed  with  values  of  E  =  10.  3  x  10®,  =  40.  5  x  10®, 

and  n  =  7.  A  solution  by  another  method  is  given  in  Reference  6.  A 
comparison  of  the  stresses  obtained  by  the  two  methods  is  shown  in 
the  table  in  Figure  2.  The  tip  displacement  is  shown  as  a  function  of 
the  load  in  Figure  3.  The  displacements  obtained  by  the  two  methods 
agree  so  closely  that  both  solutions  are  represented  by  a  single  curve. 

2.  SHEAR  LAG  SPECIMEN 

As  a  second  example,  solutions  are  obtained  for  the  shear  lag 
specimen  tested  in  Reference  3  shown  in  Figure  4.  Two  solutions  are 
obtained.  In  the  first  solution  the  same  idealization  of  the  structure 
is  used  as  that  used  in  Reference  3  (see  Figure  5).  The  second 
solution  is  found  using  the  idealization  shown  in  Figure  6.  Values  of 
the  Ramberg-Osgood  constants  of  E  -  10.  2  x  10®,  <ji  =  46.  6  x  10®, 
and  n  =  10  were  used.  These  correspond  to  the  values  for  the  R02 
stress-strain  curve  of  Reference  3. 
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ELEMENT 

FORCE,  in  kips 

Ref.  6 

Present  Analysis 

1 

11.  26 

11.  18 

2 

6.  23 

6.  19 

3 

-  6.  27 

-  6.  31 

4 

-11. 24 

-11.  11 

5 

-  0.  51 

-  0.  52 

6 

3.  35 

3.  33 

7 

6.  91 

6.  95 

8 

-  5.  59 

-  5.  55 

9 

-  4.  15 

-  4.  17 

10 

4.  47 

4.  44 

Figure  2.  Nonlinear  Truss 
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x 


Displacement  X  in  the  Nonlinear  Truss 


Figure  3.  Load  Versus  Tip  Displacement  -  Nonlinear  Truss 
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Figure  5.  MIT  Shear  Lag  Specimen  -  Configuration  I 


Figure  6.  MIT  Shear  Lag  Specimen  -  Configuration  II 
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Load  —  lb  x  I03 

Figure  7.  Comparison  of  Test  Results  of  Shear  Lag 
Specimen  with  Finite  Element  Analyses, 
Axial  Strain  at  Center  of  Stiffener 
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A  comparison  of  the  results  of  these  solutions  with  the  test  data  for 
axial  strain  at  the  center  of  the  stiffener  is  given  in  Figure  7.  While 
the  second  solution  is  in  better  agreement  with  the  test  results  than  the 
first,  the  agreement  at  large  values  of  the  load  is  poor.  A  still  more 
refined  idealization  would  probably  improve  the  solution  but,  as  shown 
in  Reference  5,  much  of  the  discrepancy  is  due  to  the  inadequacy  of  the 
tensile  stress -strain  data  in  the  range  of  large  strains. 

3.  PERFORATED  STRIP 

Theocaris  and  Marketos  (Reference  8)  obtained  results  for  a 
linear  strain-hardening  aluminum  strip  with  a  ratio  of  hole  diameter 
to  strip  width  of  1:2.  The  material  properties  were 

yield  stress  =  24.  3  kg/mm^ 

plastic  modulus  =  225.  0  kg/mm^ 

Young's  modulus  =  7000.  0  kg/mm“ 

The  finite  element  idealization  of  the  test  specimen  is  shown  in 
Figure  8;  116  nodes  and  172  triangular  elements  were  used.  A 
comparison  of  measured  and  computed  values  of  the  maximum  strain 
in  the  y  direction  at  the  edge  of  the  hole  is  given  in  Figure  9.  The 
agreement  between  theory  and  experiment  is  fairly  good.  The  same 
test  is  used  by  Marcal  and  King  (Reference  7)  for  comparison  with 
the  results  of  their  analysis  and  about  the  same  degree  of  agreement 
is  obtained. 
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Figure  8.  Perforated  Strip  Finite  Element  Idealization 
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Figure  9.  Maximum  Values  of  E  e  /crQ  for  Perforated  Strip 
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SECTION  IV 

DESCRIPTION  OF  COMPUTER  PROGRAM 

1.  INTRODUCTION 

The  FORTRAN  IV  program  for  the  elastic -plastic  analysis  of 
plane  structures  composed  of  bar  and  triangular  plate  elements  is 
described  in  this  section.  The  correspondence  between  the  program 
variables  and  the  stress-strain  law  parameters  for  each  of  the  three 
laws  available  is  given  in  Table  I. 


TABLE  I 

CORRESPONDENCE  BETWEEN  PROGRAM  VARIABLES 
AND  STRESS-STRAIN  LAW  PARAMETERS 


Stress-Strain  Law 

Program  Variables 

ILAW 

E 

EE1 

EE2 

PRR 

Ramberg-Osgood 

1 

E 

al 

n 

V 

Goldberg-Richard 

2 

E 

au 

n 

V 

Bilinear 

3 

E 

ay 

E 

1 

V 

2.  INPUT  DATA  AND  DESCRIPTION  OF  OUTPUT 

The  geometry  of  the  structure  is  determined  by  specifying  the  x 
and  y  coordinates  of  each  node  with  respect  to  a  fixed  set  of  coordinate 
axes  and  the  thickness  (cross-sectional  area  in  the  case  of  bars)  of  the 
elements.  Up  to  225  nodes  and  400  elements  can  be  handled.  The 
program  uses  a  subroutine  for  the  solution  of  simultaneous  equations 
in  band  form  written  by  Professor  E.  L.  Wilson  of  the  University  of 
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California.  Great  economies  in  storage  requirements  and  in  time 
required  for  solution  are  achieved  in  this  way;  however,  the  bandwidth 
of  the  equations  defined  by  the  idealization  of  the  structure  is  limited  in 
size.  To  meet  this  limitation  on  bandwidth  the  difference  between  the 
node  numbers  on  any  element  must  be  9  or  less.  Instructions  for  in¬ 
creasing  the  bandwidth  are  given  in  Table  II. 

The  displacement  components  in  the  x  and  y  direction  can  be 
specified  at  any  node  or  a  node  can  be  required  to  move  along  a  line 
with  a  specified  slope.  Boundary  conditions  can  be  specified  at  up  to 
29  nodes. 

The  x  and  y  components  of  load  can  be  specified  at  any  node. 
Distributed  loads  must  be  treated  as  concentrated  at  the  nodes. 

The  number  of  equal  increments  (steps)  into  which  the  applied 
loads  and  specified  displacements  are  to  be  divided  is  specified  as 
input.  It  is  also  necessary  to  specify  the  number  of  the  increment 
at  which  the  solution  is  to  start.  For  example,  if  a  number  of  incre¬ 
ments  NDIV  =  20  is  specified  and  a  value  of  the  starting  increments 
KSTART  =  5  is  used,  one  quarter  of  the  load  (displacement)  will  be 
applied  in  the  first  step,  the  rest  in  15  equal  increments.  If  it  is 
desired  to  stop  the  solution  at  an  intermediate  step  a  value  of  KSTOP 
may  be  specified.  If  the  unloading  solution  is  desired  the  value 
IUNLD  =  1  is  used. 

An  error  tolerance  must  be  specified  as  input.  After  each  cycle 
of  iteration  the  maximum  error  among  all  the  elements  is  compared 
with  the  specified  tolerance.  If  the  tolerance  is  met  the  next  load 
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TABLE  II 

PROGRAM  MODIFICATION 


The  bandwidth  is  governed  by  the  difference  between  the  node  numbers 
of  a  given  element.  The  largest  such  difference  J  determines  the 
bandwidth  in  this  program  by  the  formula  MBAND  =  2  *  J  +  3.  This 
number  cannot  be  greater  than  22  in  the  present  program.  This  is  a 
rather  small  bandwidth,  but  it  allows  225  nodes  in  a  32K  core  machine. 
To  change  the  applicable  problem  size  of  the  program  in  terms  of  the 


basic  problem  size  parameters  the 
changed: 

To  Change: 

1.  Number  of  materials 

2.  Number  of  elements 

3.  Number  of  nodes 

4.  Bandwidth 

5.  Number  of  nodes  with 

boundary  conditions 


following  dimensions  have  to  be 


Change  Dimensions  of: 

EE,  EE1,  EE  2,  PRR,  and  TAB  in 
common  statements  (presently  10) 

II,  12,  13,  14,  NTYPE,  Z,  SEF, 
SET,  EEP,  EXPL,  EYP,  EXYP  in 
common  statement,  JX  in  main 
program,  and  modify  equivalence 
statements  containing  JX 
(presently  400) 

B,  X,  XCORD,  Y,  ICODE,  FP,  F 
in  common  statements,  FE  in  main 
program,  ICODE  in  subroutine 
DCODE,  XX  in  subroutine  ELEM 
(presently  225  or  2  x  225  =  450) 

B  in  common  statements 
(presently  22) 

BC  in  common  statements 
(presently  30) 
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increment  is  applied,  if  not,  the  iteration  is  continued.  If  the  tolerance 
on  error  is  not  met  when  the  allowable  number  of  iterations  is  reached 
the  solution  is  stopped.  A  detailed  description  of  the  input  data  format 
is  given  in  Table  M. 

TABLE  III 

INPUT  DATA  FORMAT* 


Card  1  TITLE  CARD  (72H) 

Col  1-72  Any  alphanumeric  information 

Card  2  PROPERTIES  CARD  (1415) 


Col 

1-  5 

NNODE 

- 

number  of  nodes  (maximum  225) 

6-10 

NELEM 

- 

number  of  elements  (maximum  400) 

11-15 

ILAW 

- 

1  Ramberg-Osgood  Law 

- 

2  Goldberg -Richard  Law 

- 

3  Bilinear  Law 

16-20 

IUNLD 

- 

1  Unloading  following  loading 

- 

0  Loading  only 

21-25 

MAT 

- 

number  of  materials  used  (maximum  10) 

26-30 

MAXBND 

maximum  bandwidth,  MAXBND  =  22 
for  this  program 

31-35 

NBC 

number  of  boundary  conditions  with 
prescribed  displacement.  The 
maximum  number  is  30  in  this 
program. 

Card  3 

MATERIAL  PROPERTIES  CARDS  (E15.  8,  3F10.  5) 

Col 

1-15 

EE 

- 

modulus  of  elasticity 

16-25 

EE1 

secant  yield  stress,  ultimate  stress, 
yield  stress 

26-35 

PRR 

- 

Poisson's  ratio 

36-45 

EE  2 

shape  parameter,  plastic  modulus 

Card  4 

CONTROL  CARD 

(615 

,  F10.  0) 

Col 

1-  5 

NDIV 

- 

number  of  load  increments 

6-10 

NIT 

- 

maximum  number  of  iterations  per  step 

11-15 

NPRINT 

print  output  for  each  NPRINT  incre¬ 
ment.  (e.  g. ,  if  NPRINT  =  3,  for 
increments  3,  6,  9  etc. ) 

16-20 

KSTART 

— 

number  of  increments  at  which 
solution  is  to  start 

21-25 

KSTOP 

— 

number  of  increments  at  which 
solution  is  to  stop 

26-30 

NLOAD 

“■ 

number  of  nodes  at  which  loads  are 
specified 

31-40 

TOL 

- 

error  tolerance 
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TABLE  III  (CONTD) 


Card 

5 

NODE 

CARDS  (415. 

,  5F10.  0) 

Col 

1-  5 

Node  number 

6-10 

IBCX  = 

1, 

if  displacement  in  x -direction  is 
specified 

11-15 

IBCY  = 

1, 

if  displacement  in  y-direction  is 
specified 

16-20 

IBCS  - 

1, 

if  slope  is  specified 

21-30 

XCORD 

- 

x  coordinate  of  the  node 

31-40 

YCORD 

- 

y  coordinate  of  this  node 

41-50 

BCI 

- 

specified  displacement  in  x-direction 

51-60 

BC2 

- 

specified  displacement  in  y-direction 

61-70 

BC3 

- 

specified  slope  at  the  node 

Card 

6 

ELEMENT  CARDS 

(515,  F10.  0) 

Col 

1-  5 

Element  number 

6-10 

11 

- 

nodes  defining  the  element 

11-15 

12 

- 

nodes  defining  the  element 

16-20 

13 

- 

nodes  defining  the  element 

21-25 

NTYPE 

- 

material  type 

26-35 

Z 

- 

element  thickness  or  cross-sectional 

area 

Card 

7 

LOAD 

CARDS  (15, 

2F10.  0) 

1-  5 

Node  number 

6-15 

x-component  of  force 

16-25 

y-component  of  force 

♦NOTE:  Input  data  information  in  Table  III  is  self  explanatory.  The 
use  of  more  than  one  material,  however,  may  need  some  clarification. 
The  number  of  materials  "MAT"  specified  in  the  field  of  card  2  defines 
the  number  of  material  properties  cards.  The  sequencing  of  these 
cards  in  turn  defines  "NTYPE"  in  the  element  card,  for  example,  if 
the  element  uses  material  specified  in  the  second  material  properties 
card,  integer  2  is  placed  in  the  field  corresponding  to  "NTYPE". 


27 


AFFDL-TR-68-39 


The  nodal  forces  and  displacements,  the  maximum  error  and  the 
number  of  the  element  in  which  it  occurs  are  printed  out  at  the  end  of 
each  step  (increment).  The  cartesian  components,  principal  values, 
and  direction  of  stress  and  strain  are  printed  out  at  the  user's  option 
by  specifying  a  value  of  NPRNT  as  input.  For  example,  a  value  of 
NPRNT  =  3  will  cause  the  stresses  and  strains  to  be  printed  out  for 
increment  numbers  divisible  by  three.  The  directions  of  the  principal 
axes  of  stress  are  defined  by 


- 


± 

2 


ton 


4> 


< 


7T 

2 


The  value  of  (f>  in  degrees  is  printed  out.  In  the  case  of  strain  the 
principal  directions  are  defined  by 


4>--  ~ 


2 


ton 


7 r 
2 


This  value  is  also  printed  out  since  in  general  the  principal  axes  of 
stress  and  total  strain  do  not  coincide  when  plastic  flow  has  taken 
place. 


The  effective  stress  and  the  effective  plastic  strain  are  also  given 
as  output  for  each  element. 

The  input  data  is  printed  out  at  the  start  of  the  program  to  aid  in 
problem  identification  and  checking. 
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APPENDIX 

COMPUTER  PROGRAM  LISTING 


The  FORTRAN  IV  Source  Program  and  three  sample  data  cases  are 
listed.  The  first  case  is  the  nonlinear  truss  with  Ramberg-Osgood 
representation  of  one  material.  This  data  is  associated  with  Figures  2 
and  3.  The  second  case  is  the  same  truss  problem  slightly  changed  to 
show  the  introduction  of  more  than  one  material.  The  third  case  is 
Configuration  I  for  the  MIT  test  specimen. 

The  source  deck  of  the  computer  program  described  herein  can  be 
obtained  by  contacting  AFFDL  (FDTR/BERKE),  WPAFB-Ohio,  45433. 
(513-25-53418). 
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SIBJOB 

SIBFTC 

C 


c 

c 

c 


VER7 

ELASTIC  PLASTIC  FINITE  ELEMENT  PROGAM 
WITH  THREE  STRESS  STRAIN  LAw  OPTIONS 
COMMON /AOD/  Efc ( 10 ) f E E 1 ( 10 ) , EE 2 ( 10 ) ,  PRR  (  10) 

COMMON  EfCCtGlfE2,PR.EPR,X21fY21fX31tY31fX32fY32fXERRf 

1  N2 , NE  L  EM , K  E  L , I L A W , MAT , N BC , 

2  B( 450, 22 ) ,BC ( 30,3) , TAB( 101, 20) , IFl X( 2) , 

3  X ( 4  50 ) v  XCORO (  225) ,Y(225)»ICODE(225  ), 

4  FP(450)  ,F(450)  » 

5  I  1 (400 )  ,  i 2 i 400) f  13 ( 400) 9NTYPE (400)  9  7(400) 9  14(400)  , 

6  SEF(40G),SET (400) ,EEP( 400) tEXPL(400) fEYP(400) ,EXYP 

7  NNOOE  f  MB AND 
DIMENSION  JX(400,3)  ,FE(450) 

EQUIVALENCE  (  J  X  9  ID  9  (  JX  (  40 1  )  9  I  2  )  9  (  JX  (  801)  9  13  ) 

EQUIVALENCE  ( IFIX( 1) v ieCX) , ( IFIX(2)f I8CY) 

****  READ  AND  PRINT  DATA  **** 

10  READ  ( 5 9  20 ) 

I  ONE  =  1 

20  FORMAT ( 72H  BCO  INFORMATION 
1  ) 

WRITE  (6,30) 

30  FOKMATt 1H1 ) 

WRITE  (6,20) 

RE  AO ( 5,40 )  NNOOE , NE LEM , I  LAW, I UNL 0 , MA T , MA XBN D , NBC 
40  FORMAT ( 1415) 

READ( 5, 50)  ( EE (  I  )  , EE  1  (  I  )  ,PRR( I ) , EE 2 ( I ) ,  I  =  1 , MAT ) 

50  FORMAT (E15«8,3F10.5) 

RE AD( 5, 60 ) NDl V,NIT ,NPRNT , K START , KS TOP , NLOAD , T OL 
60  FORMAT ( 6I5,F10«G) 

I F ( K STOP. EC ,0)  K  STOP=NO  l  V 
GO  T0(70,90, 110) , ILAW 

70  MR  I TE( 6,80)  (  I  ,  EE ( I  ) , EE  1 ( I ) , EE2 ( I) ,PRR (I ) , TOL , I =1 ,MAT ) 

80  FORMAT! 1H014X 1 8HRAMBERG  OSGOOD  LAW/ 

115X,  30HMATER  IAL -  13/ 

215X30HMO0ULUS  OF  ELASTICITY - E12.4/ 

3 1 3X30HS  EC  ANT  YIELO  STRESS - E12.4/ 

4  1  5X  3  OHS  HA  PE  PARAMETER -  E12.4/ 

515X30HP0ISS0N*  S  RATIO - F8.4/ 

6 1  5X3CHE RROR  TOLERANCE - F8.4) 

GO  TO  130 

90  WRITE(6,  100)  (  I,EE(I),EE1C  I),EE2(I)  ,PRR(  I  )  , TOL, 1  =  1 , MAT ) 

100  FORMAT ( 1H014X20HGOL06ERG  RICHARO  LAW/ 

115X,30HMATERIAL -  13/ 

215X30 H MODULUS  OF  ELASTICITY - E12.4/ 

315X30HULT  IMATE  STRESS - E12.4/ 

415X30HSHAPE  PARAMETER - E  1  2 . 4/ 

515X30HP0I  SSGNT  S  RATIO - F8.4/ 

61 5X30HERR0R  TOLERANCE - F8.4) 

GO  TO  130 

110  WRITE (6, 120)  (  I,EE(  I  )  ,EE1(I),EE2(I) ,PRR  l  I  ) , TOL ,1  =  1 9  MAT ) 

120  FORMAT! 1 HO  14X12 HBILINEAR  LAW/ 

1 1  5X,  3CHMAT  ERIAL -  13/ 

2 1 5X3CHM0DULUS  OF  ELASTICITY - E12.4/ 

315X30HY  IELO  STRESS - E  12 .4/ 

415X30HPLASTIC  MODULUS - E  12 . 4/ 

515X30HPOISSON'  S  RATIO - F8.4/ 

615X30HE  RROR  TOLERANCE - F8.4) 

130  WRITE  ( 6  9 140 ) NNOOE  9NELEM9NOIV9NIT 

140  FORMAT ( 1H014X30HN0.  OF  NODES  NNOOE  =14/ 1 5X30HNO. 


VER70000 

V  ER7000 1 
VE R  70002 
VER70003 
VER70004 

V  ER70005 
VER70006 
VER70007 
VER70008 
VER70009 

(400),  VER70010 

VER700 1 1 
VER70012 
VER70013 
VER700 14 
VER700  1  5 
VER700 16 
VER70017 
VER70018 
VER700 1 9 
VER70020 
VER7002 1 
VER70022 
VER70023 
VER70024 
VER70025 
VER70026 
VER7002  7 
VEK70028 
VER70029 
VER70030 
VEK70031 
VERY 0032 
VER70033 
VFR700  34 
VER70035 
VER70036 
VER70037 
VER70036 
VER  70039 
VER70040 

V  ER  70  04 1 
VFR70042 

V  ER  7  0  04  3 
VER70U44 

V  ER70G45 
VEK70046 
VER70047 
VER70048 

V  ER  7004  9 
VER  70050 
VEK70051 
VER70052 

V  E  K  700  53 
VEK70054 
VER70055 
VER70056 
VER70057 
VER70058 
VER70059 

OF  E  L  EM V  ER70060 
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1  ENT S  NELEM  =1 4/ 1 5X30HN0 .  OF  STEPS 

2HN0 »  OF  ITERATIONS/STEP  NIT  =14) 

DO  150  1=1, NBC 
DO  150  J=  1 , 3 
150  BC( It J)=0 
IC=1 

WRITE (6, 160) 

160  FORMAT ( 25H0B0UNDARY  CONDITION  ARRAY/ 10H0  NODAL  PT 1 5X 1HX23X 1HY 
120X7HSLIDING/1H  14X4HC0DE7X  5H VALUE 9X4HC0DE  7X5HV ALUE9X4HC ODE 
27X5HV ALUE ) 

****  NODE  COORDINATES  AND  BOUNDARY  CONDITIONS  **** 

DO  200  J  =  1 , NNODE 

READ(5, 170)  K ,  I 8C X , I 8C Y , I 3CS , XCORD < K ) ,Y(K) ,BC1,8C2,BC3 
170  FORMAT!4I5,5F10.0) 

I F  (  IBCX  +  IBCYMBCS.NE.O)  WRI TE < 6, 180) K , I BCX , 8C1 , I BCY, BC2 ,  IBCS,8C3 
180  FORMAT! I7,3X,3! I8,1PE17.7) ) 

ICGDE(K)=IBCS+10*IBCY+100*IBCX 

IFIBC1+BC2+BC3.EQ.0. )  GO  TO  200 

IC0DE!K)=IC0DE!K)+IC*1000 

BC (  IC,  1  )=BC1 

BC(  I C  *  2  )  =  BC2 

BC (  IC, 3)=BC3 

IC=IC+1 

I F ( IC.LE.NBC )GO  TO  200 
WRITE16, 190) 

190  FORMAT ( 54H0  MORE  THAN  29  NODES  HAVE  NON  ZERO  BOUNDARY  CONDITIONS) 
GO  TU  620 
200  CONTINUE 


****  ELEMENT  PROPERTIES  **** 


210 


READ! 5,210) (K, I1(K),I2(K),I3(K) , NTYPE ( K ) 
FORMAT!  5I5,F10.0) 


Z (K)  , J  =  l, NELEM) 


**** 


LOADS  **** 


220 


230 

240 

250 

260 


• 1 ) , F ( 2* J  ) 


N2=2*NN0DE 
DO  220  K= 1 , N2 
F  (  K ) =0 

I F ( NLOAD. EG • 0 )  GO  TC  250 
DO  230  K= 1 , NLOAD 
READ! 5, 240)J ,F (2*J 
FORMAT!  I  5, 2F 10.0 ) 

CONTINUE 

WRITE! 6, 260) IK,XCURDIK) ,F(2*K-1) ,  Y  !  K )  , F { 2*K ) , I  CODE ! K ), K  =  1 , NNODE ) 
FORMAT! 10H0  NODAL  PT8X7HX-C00RD8X7HX-F0RCE8X7HY-C00RD8X7HY-F0RCE 
111X4HC0DE//(4X, I3,5X4F15.4, 115) ) 

WRI f  b ( 6 , 2  70 ) 

NODE  1  NODE  2  NODE  3 
MATERIAL  TYPE,//) 


270  FORMAT! 1H0/// 1GX ,90HELEMENT 
IT  TYPE  AREA  OK  THICK. 

DO  300  K=l, NELEM 

IF! I  3 ! K )  .EQ.O)  WRITE!6,2  80)  K,I1(K),I2(K),I3(K),Z(K)  , NTYPE ! K ) 
IFII31KJ.NE.G)  WR I T  E ( 6 , 2  90 )  K,I1(K),I2(K),I3!K),Z(K),NTYPE!K) 
FORMAT (  6X,4I9, 11X,3H8AR, E22.5, 114) 

6X,4I9, 11X,5HPLATE,E20.5, I  14) 


280 

290 

300 


FORMAT ( 
CONTINUE 


INITIALIZATION 


310  XD I V  =ND I V 


NDIV  =I4/15X30VER70061 
VER70062 
VER70063 
VER70064 

V  ER70065 
VER70066 
VER70067 
VER70068 
VEK70069 
VER70070 
VER70071 
VER70072 
VER70073 
VER70074 
VER7007  5 
VER70076 
VER70077 

V  ER70078 
VER70079 

V  ER70080 

V  ER70081 
VER70082 

V  ER70083 
VER70084 

V  ER70085 
VER70086 
VER70087 
VER70088 

V  ER70089 
VER70090 
VER70091 
VER70092 
VER70093 

V  ER70094 

V  ER70095 
VER70096 
VER70097 
VER70098 
VER70099 
VER70100 
VER70101 
VER70102 
VER70103 
VER70104 
VER70105 
VER70106 
VER70107 
VER70108 
VER70109 
VER70110 

ELEMEN V  ER70 1 1 1 
VER70112 
VER701 13 
VER70114 
VER70115 
VER70116 
VER70117 
VER70118 
VER70119 
VER70120 
VER70121 
VER701 22 
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GO  T0(320f  340t  340) , I  LAW 
320  CONTINUE 

00  330  1=1, MAT 
E  =  EE (  I  ) 

E  1  =  E  E  1  (  I  ) 

E  2=E  E2 ( I ) 

CC=E 1 / E 

Gl=(7.0*E/3.0)**( 1.0/E2)*E1**( 1.0-1. 0/E2) 
CALL  TABLE(I) 

330  CONTINUE 
340  CONTINUE 

****  OETERMINE  BANO  WIDTH  **** 


350 


360 

370 

380 


390 


400 


410 


00  350  K= 1 , NELEM 
I  4  (  K  )  =  I  3  (  K  ) 

I F ( 13 ( K ) .EO. 0)  JX(K,3)=JX(K, 1 ) 

J  =  0 

00  380  N=l, NELEM 
00  380  1=1,3 
DO  370  L= 1 , 3 

KK= I ABS ( JX ( N , I ) - JX ( N , L ) ) 

I F ( KK- J )370, 370, 360 
J  =  KK 

CONT INUE 
CONTINUE 
MB AND=  2* J ♦ 3 
I F ( MBAND.GT • MAX  BND ) 

IF(MBANO«Gr  *MAXBN0) 

FORMAT ( 1H010X20HBAN0 
DO  400  1=1, NELEM 
I  3  (  I  )  =  I  4(  I  ) 

00  410  1=1, N2 
00  410  J= 1 , MBANO 
B (  I  ,  J  )  =0 • 

CALCULATION  OF  STIFFNESS  MATRIX 
CALL  STIFF 


WR I TE ( 6 »  390 )  MBAND 
GO  TO  10 

WIOTH  TOO  L ARGE5X  6HMB AND= I  4 ) 


****  REDUCE  MATRIX  **** 

CALL  SYMSOL ( 1 > 

****  INCREMENT  LOAOS , AOD  PLASTIC  FORCES  ANO  SOLVE  FOR 

DO  420  1=1, NELEM 
SET (  I  ) =0 
SEF ( I ) =0 
EE P (  I  >  =0 
EXPL ( I )=C 
EY P (  I  )  =0 
420  EXYP ( I ) =0 
KO=KST  ART- 1 
KU  =  KO 

DO  430  1=1, N2 
FP (  I  ) =0 
GO  TO  490 


430 

440 

450 


VER  70 1 23 
VER70124 
VER70125 
VER70126 
VEK70127 
VER70128 
VER70129 
VER70130 
VER70131 
VER70132 
VER70133 
VER70134 
VER70135 
VER70136 
VEK70137 
VER70138 
VER  70 1 39 
VER70140 
VER70141 
VER70142 
VEk70 143 
VER70144 
VEK70145 
VEK70146 
VER70147 
VER70148 
VEK70149 
VER70150 
VER70151 
VER70152 
VER  701 53 
VER 7 0154 
VER70155 
VER70156 
VER70157 
VER70158 
VER70159 
VER70160 
VER701 61 
VER70162 
VER70163 
VEK70164 
VER70165 
VER70 166 
DI$PLACEMENT$**VER70167 
VER70168 
VFR70169 
VERIO  170 
VER70171 
VEK70172 
VER70173 
VER70174 
VEK70175 
VER70176 
VER70177 
VEK70178 
VEK70179 
VEK70180 


WRITE  ( 6,450)KU» ( I ,FE( 2*1  —  1 ) , F  E ( 2* I ) , X ( 2* I - 1 ) , X ( 2* I ) , I  =  1 , NNODE )  VER70181 

FORMAT ( 1H120X38HFORCES  AND  DISPLACEMENTS  FOR  INCREMENT , I4//10X4HN0VER70182 
1 DE  5X  7HX-F0RCE8X7HY- FORCE  9X  8HX-0ISPL.7X8HY-DISPL./(9XI3,2F15.3,5X2EVER70183 
215.4  ))  VER70184 
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WRITE  ( 6,460)XhRR,KEL, IT 

460  FORMAT! 13H0MAX.  ERROR  = F 8 . 5 , 5X 14H I N  ELEMENT  NO. I  4 , 5X 17HN0. 
1  AT  I  ON S I  4 ) 

470  I F ( MOD ( KO »  NPRNT ) 1490,4  80,490 
480  CALL  OUTPT 

I  F  (  KU  •  F(v  •  K  S  TOP  1  GO  TO  620 
GO  TO  600 

490  IF(KO.EQ.KSTOP)  CALL  OUTPT 
IF(K0.EG.K$TUP)G0  TO  620 
500  K0=KG+  I  ONE 
KU=KU+ 1 

IF( KG-NDIV) 510, 510,620 
510  XK0  =  K 0 

DO  520  1=1, N2 
520  FE (  I )=XKO/XDIV*F( I  ) 

DO  530  K= 1 , NELEM 
530  SEF(K)=$ET(K) 

I  T  =  0 

540  DO  570  1=1, NNODE 

I F ( I  COD  E ( I )  •  EQ  •  0 )  GO  TO  570 

CALL  DCOPF (  I  CODE, I , IBCS, IBCX,  IBCY, IC, IX,  IY,NBC) 

I F ( IBCS.NE.l )  GO  TO  550 
ALF  =  BC (  I  C , 3  ) 

F  P (  I  X  )  =  F  P  (  IX  ) + A  L  F*  F  P (  IY) 

FP( I Y ) =0. 

550  DO  560  N= 1 , 2 

I F ( IFIX(N) .NE.l)  GO  TO  560 
I R= IX+N-l 
FP( IK)=0. 

560  CONTINUE 
570  CONTINUE 

****  SOLVE  FOR  DISPLACEMENTS  **** 

DO  580  1=1, N2 
580  X(  I  )  =  E E (  I ) +E P (  I  ) 

CALL  SYMSOL ( 2 ) 

CALCULATE  TOTAL  STR A  I  NS , STRESSL S  AND  PLASTIC 
FORCES  AND  STRAINS  FUR  EACH  ELEMENT 

DO  590  1=1, N2 
590  FP(  I  1=0 
XERR=0 • 0 
KEL  =  U 

CALL  STRAIN 

****  PICK  LARGEST  ERROR  AND  DETERMINE  wHEN  TO  REITERATE  **** 

I  T  =  I  T  +  1 

I F ( X  E  R  R - T  0  L 1440,440,600 
toOO  IF( I T-NIT 1540,610, 610 
610  K  0  =  N  b  IV 

iunld=o 

GO  TO  440 

620  IF( IUNLD.EQ.O)  GC  TO  10 
I  U  N  L  D  =  0 
IUNE=-1 
K  S  T  (J  P  =  0 
GO  TO  44C 
END 

SIGFTC  ELM 


VER70185 
OF  I  TER VER70 1 86 
VER70187 
VER 70188 
VER70189 
VER70190 
VER70191 
VER70192 
VER70193 
VER70194 
VER701 95 
VER70196 
VER70197 
VER70198 
VER70199 
VER70200 
VER70201 
VER70202 
VER70203 
VER70204 
VER70205 
VER70206 
VER70207 
VER70208 
VER70209 
VER70210 
VER70211 
VER70212 
VER702 1 3 
VER70214 
VER702 1 5 
VER70216 
V6R70217 
VER70218 
VER70219 
VER70220 
VER70221 
VER70222 
VER70223 
VFR702  24 
VER70225 
VER70226 
VER70227 
VER70228 
VER70229 
VER70230 
VER7023 1 
VLR70232 
VER70233 

V  ER70234 
VER7023  5 
VEK70236 
VER  702  37 
VER  702  3  8 
VER70239 

V  ER70240 
VE R  7024 1 

V  ER70242 
VER70243 

V  ER70244 
VER70245 
ELM  0000 
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SUBROUTINE  ELEM \ M ) 

COMMON/ ADD/  EE(10),EE1(10),EE2(10) ,PRR( 10) 

COMMON  EfCC»GlfE2*PRfEPRfX21fY21fX31fY31fX32fY32fXERRf 

1  N2 , NEL  EM  ,  KE  L ,  I L AW , M AT , NBC , 

2  B(450,22),BC(30,3) ,TAB( 101,20)  ,IFIX(2)  , 

3  X( 450) , XCORD( 225) ,Y ( 225 ) , I CODE ( 225 ) , 

4  FP(450) *F(450) , 

5  11 (400 )  ,  12(400 ) , I 3( 400) ,NTYPE(400)  ,Z(400) *  14(400) * 

6  SEF(400) , SET (400  WEEP (40G) ,EXPL( 400) ,EYP( 400) ,EXYP( 400) 

7  NNODE «  MbAND 
DIMENSION  X  X (  22  5 ) 

EQUIVALENCE  <XCOKD,XX) 

J 1= I  1 ( M  ) 

J  2= I  2 ( M  ) 

J  3=  I  3  (  M  ) 

X21=XX( J2)-XX( Jl) 

Y21=Y( J2)-Y< Jl) 

IFU3.EQ.0)  GO  TO  10 
Y32  =  Y ( J  3 ) - Y ( J2 ) 

Y  3  1  =  Y  (  J  3  )  -  Y  (  J  1  ) 

X32-XX(J3)-XX(J2) 

X31=XX( J3)-XX( Jl) 

RETURN 

10  Y32=SQRT( X2 1**2+Y2 1**2 ) 

RETURN 

END 

SIBFTC  DCUD 
CDCOD 

SUBRUUT INE  DCODE( I  CODE* I  * IBCS* I BCX *  I BC Y * IC *  I  X , I Y, NBC ) 

DIMENSION  ICODE ( 225 ) 

I  BC  S  =  MOD (  I  CODE (  I ) , 10) 

I  BCX  =  MOD ( I  CODE ( I ) *  1000) /IOC 
I  BC Y  =  MOD ( I  CODE (  I ) *  100)/ 10 
IC=MOD(  ICODE ( I ) , 100000)/ 1000 
I  X  =  2* 1-1 
I  Y  =  I  X  +  1 


I  F (  IC.EQ.O)  I C  =  NBC 

RETURN 

END 

SIBFTC  PNEW 

C2222  PLASTIC  STRAIN  DETERMINATION 

SUBROUTINE  PNEwl  ( E EPK * EET * E * E 1 , E2 ) 
J  =  1 

XU=EE  T 
XL  =  0 

10  EEPK=*5*( XL+XU ) 

20  Y=EEPK+E1/E*EEPK** ( 1.0/E 2) -EET 
3  0  YP  =  1 .0+El/E/E2*EEPK**( 1.C/E2-1 .0) 

J  =  J  +  1 

I  F ( J-50 ) 40  *  40  *  100 
40  IF (Y ) 50, 100,60 
50  XL=EtPK 
GO  TO  70 
60  XU=  E EPK 
70  XT=EtPK-Y/YP 

IF (XU-XT ) 10, 10,80 
80  IF ( XT-XL ) 10, 10, 90 
90  EE PK  =  X T 

DIFF=ABS(Y/YP/EEPK) 

I F ( D I FF- *00001 ) 100,100,20 
100  RETURN 
END 


ELM  0001 
ELM  0002 
ELM  0003 
ELM  0004 
ELM  0005 
ELM  0006 
ELM  0007 
ELM  0008 
ELM  0009 
ELM  0010 
ELM  0011 
ELM  0012 
ELM  0013 
ELM  0014 
ELM  0015 
ELM  0016 
ELM  0017 
ELM  0018 
ELM  0019 
ELM  0020 
ELM  0021 
ELM  0022 
ELM  0023 
ELM  0024 
ELM  0C25 
ELM  0026 
DCODOOOO 
DCUDOOO 1 
DC0D00G2 
DCOD0003 
DCOD0004 
DCGD0005 
DCCU0006 
DCQU0007 
DCODOO  08 
DC000009 
DCODOO 1 0 
DC0D0011 
DCODOO 1 2 
P  NEWOOOO 
PNEW0001 
PNEW0002 
PNEW0003 
PNEW0004 
PNEW0005 
PNEW0006 
PNE WOOD  7 
PNEH0GC8 
PNfcWO  009 
PNEWOOIO 
PNE  WOO  1 1 
P  NtWOO 1 2 
PNEWOO 1 3 
PNh WOO  1 4 
P  NE  WOO  1 5 
PNEWOO 1 6 
PNEWOO 1 7 
PNEWOO 1 8 
PNEWOO 1 9 
PNEW0020 
PNEW0C21 
PNEW002  2 
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SIBFTC  PLSTR 

CTABL  STRAIN-PLASTIC  STRAIN  TABLE 
SUBROUTINE  TABLE ( K ) 

COMMDN/ ADD/  EE(IO),EEI(IO) ,EE2( 10) ,PRR(10) 

CDMMDN  E,CC,G1,E2,PR,EPR,X21, Y21,X31 ,Y31,X32,Y32,XERR, 

1  N2,NELEM,KEL, I  LA W, MAT ,NBC, 

2  B(450,22),BC(30,3) ,TAB( 1 01 , 20 ) ,  IF  I X ( 2 )  , 

3  X ( 450 ) t  XCDRD (225),Y(225) * ICODE<2  25) t 

4  FP(450) ,F(450) , 

5  11(400) , 12(400) , 13(400) ,NTYPE(400) ,Z(400) , 14(400) t 

6  SE  F ( 400 ) , SET (400) ,EEP(4C0) , EXPL(400) , EYP(400) ,EXYP(400) , 

7  NNODE  v  MB AND 
I  1 2  =  2*K 

I  1 1=112-1 
T  AB  (  1 f 1 1 1 ) =0 . 

T AB ( It  I  1 2 ) =0  • 

DO  20  1=1,101 
I F ( 1-1)20,20,10 

10  TAB( I, I I1)=FLDAT( I-l)*CC/5. 

EET  =  T AB ( 1 ,  1 1 1 ) 

CALL  PNEWKEEPK, EET,E,G1,E2) 

T AB ( I , I  I  2 ) =EEPK 
20  CONTINUE 
RETURN 
END 

SIBFTC  STIF 

SUBROUTINE  STIFF 

CDMMDN/ ADD/  EE(10),EE1(10),EE2(10) ,PRR( 10) 

COMMON  E,CC,G1,E2,PR,EPR,X21,Y21,X31,Y31,X32,Y32,XERR, 

1  N2,NELEM,KEL,ILAW,MAT, NBC, 

2  B ( 450, 22) ,BC(30, 3) ,TAB( 1 01 , 20 ) , I F I X ( 2 ) , 

3  X(450) , XCDRD (225) ,Y(225) , ICDDE (225) , 

4  FP(450) ,F(450) , 

5  II (400) , 12(400) ,  13(400) ,NTYPE( 400) , Z ( 400 ) , 1 4 ( 400 ) , 

6  SEF(400) , SET (400) ,EEP(400) ,EXPL(400) , E YP ( 400 ) , EXYP ( 400 ) , 

7  NNDDE , MBAND 

DIMENSION  FFP ( 2 ) ,  INDDE ( 3 ) , LNODE ( 6 ) ,DSK(6,6) 

EQUIVALENCE  ( IFIX(l) ,IBCX) , ( I F I X ( 2 ) , I BC Y ) 

EQU I VAL  ENCE ( INODE ( 1 ) , Jl) , ( INODE ( 2) , J2) , ( INODE (3) , J3) 

DO  80  M= 1 , NE LEM 
J Is  1 1 (M ) 

J2= I  2 ( M ) 

J3* I  3 ( M ) 

CALL  ELEM(M) 

NTY=NTYPt ( M ) 

E  =  E  E ( NT Y ) 

PR=PRR ( NT Y ) 

I F ( J 3 • EQ • 0 )  GO  TO  30 
IF(NTY.LE.MAT)  GO  TD10 
GD  TO  160 
C 

C  ****  STIFFNESS  MATRIX  CALCULATIONS  FDR  TRIANGULAR  ELEMENTS  **** 

C 


10  JM A T  =6 
AE=E*Z ( M ) 

A123=X21*Y31-X31*Y21 
A1 2 3= AB S ( A  1 23  ) 

E  T  1  =  AE/  (  2  .0*A123*  (  1  .0-PR**2)  ) 
ET2=AE/(4.0*A123*( 1.0+PR)) 

DSK (1,1)=  ET1*Y32**2  +ET2*X32**2 

DSK ( 2 , 1 )=-ETl*PR*Y32*X32  -ET2*Y32*X32 

DSK ( 2 , 2  )  =  ET1*X32**2  +ET2*Y32**2 


PLSTOOOO 
PLST0001 
PLST0002 
PLST0003 
PLST0004 
PLST0005 
PLST0006 
PLST0007 
PLST0008 
PLST0009 
PLST0010 
PLST001 1 
PLST0012 
PLST0013 
PLST0014 
PLST001 5 
PLST0016 
PLST0017 
PLST001 8 
PLST001 9 
PLST0020 
PLST002 1 
PLST0022 
PLST0023 
PLST0024 
STIFOOOO 
STIF0001 
STIF0002 
STIF0003 
STIF0004 
STIF0005 
STIF0006 
STIF0007 
STIF0008 
ST  I F 000 9 
STIF0010 
STIF0011 
STIF0012 
STIF0013 
STIF0014 
STIF0015 
STIF0016 
STIF0017 
STIF0018 
STIF0019 
STIF0020 
STIF0021 
STIF0022 
STIF0023 
STIF0024 
STIF0025 
STIF0026 
STIF0027 
STIF0028 
STIF0029 
STIF0030 
STIF0031 
STIF0032 
STIF0033 
STIF0034 
STIF0035 
STIF0036 
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DSKI 3» 1 )=-ET 1*Y31*Y32 

-ET2*X32*X31 

STIF0037 

OSKC  3»2  )  =  ET1*?R*Y31*X32 

♦ET2*Y32*X31 

STIF0038 

DSKC3,3)=ET1*Y31**2 

♦ET2*X31**2 

ST1F0039 

DSK  C  A  » 1 ) =  ET1*PR*Y32*X31 

♦ET2*Y31*X32 

STIF0040 

OSKC 4,2)=-ETl*X31*X32 

-ET  2*Y31*Y32 

STIF0C41 

DSKC4,3)=-ET1*PR*Y31*X31 

— ET  2*Y31*X31 

STIF0042 

DSK ( A  »4 ) =  ET1*X31*^2 

♦ET2*Y31**2 

STIF0043 

DSK ( 5  » 1  )=  ET1*Y21*Y32 

♦ET2*X32*X21 

STIF0044 

DSK(5»2)=-ET 1*PR*Y2 1*X32 

— ET2*Y32*X21 

STIF0045 

DSKI5»3)=-ET1*Y31*Y21 

-ET  2*X31*X2 1 

STIF0046 

DSK ( 5?4 )  =  ET1*PR*Y21*X31 

♦ET2*Y31*X21 

ST1F0047 

DSK 1  5  »  5  )  =  ET 1*Y21**2 

♦  ET  2*X2 1**2 

STIF0048 

DSKI6,1)=-ET1*PR*Y32*X21 

— ET  2*Y21*X32 

ST1F0049 

DSK ( 6»  2 )=  ET 1*X32*X21 

♦ET  2*Y21*Y32 

STIF0C50 

DSK  Co*  3 )=  ET 1*PR*Y31*X21 

♦ET2*Y21*X3l 

ST1F0051 

DSK(6»4 )=-ETl*X31*X21 

-ET2*Y21*Y31 

STIF0052 

DSK|6,5)=-ET1*PR*Y21*X21 

— ET2*Y21*X21 

ST1F0053 

DSK ( 6 1 6 ) =  ET 1*X21**2 

♦ET2*Y2i**2 

STIF00S4 

DO  20  I =1 »  JMAT 

ST1F0055 

DO  20  J  =  I »  JMAT 

ST1F0056 

20  DSK IItJ)=DSKIJ»I) 

STIF0057 

GO  TO  50 

STIF0058 

STIF0059 

****  STIFFNESS  MATRIX  CALCULATIONS  FOR  BARS  **** 

STIF0060 

STIF0061 

30  ET1=Z(M)*E/Y32**3 

STIF0062 

FFP ( 1 )  =  X2 1 

ST1F0063 

FFP(2)=Y21 

ST1F0064 

DO  40  1=1,2 

STIF0065 

DO  40  J=l, 2 

STIF0066 

DSKC I,J)=ET1*FFP(I)*FFP( 

J) 

STIF0067 

DSKI 1+2 ,J ) — -DSK ( I , J ) 

STIF0068 

DSKC ItJ+2) =-DSK ( I , J ) 

ST1F0069 

40  DSKC l+Z9 J+2)=DSK< I , J) 

STIF0070 

JHAT=4 

ST1F0071 

STIF0072 

****  INCORPORATION  OF  ELEMENT  MATRICES  INTO 

STIF0073 

COMPLETE 

STIFFNESS  MATRIX  **** 

STIF0074 

STIF0075 

50  JMAT2=JMAT/2 

STIF0076 

K=0 

STIF007r 

DO  60  1=1 »  JMAT 2 

STIF0078 

DO  60  J  =  1 1 2 

STIF0079 

K=K+1 

STIF0080 

60  LNODE l K I =2* I NODE  I  I  I  — 2* J 

ST  IF 0081 

DO  80  1  =  1  * JMAT 

ST1F0082 

KI=LNODE( I ) 

STIF0083 

DO  80  J  =  1 1 JMAI 

STIF0084 

KJ=LNODE(J) 

STIF0085 

IFCKJ-KI 180,70,70 

STIF0086 

70  K=K J-K l *1 

ST1F0087 

8 ( KI t K ) =B ( K I »K  )  +  DSK ( I,J» 

STIF0088 

80  CONTINUE 

STIF0089 

ST1F0090 

****  DISPLACEMENT  BOUNDARY 

CONDITIONS  **** 

STIF0091 

STIF0092 

DO  150  1=1 tNNCDE 

STIF0093 

IF  1 I CODEC  I ) .EO.C)  GO  TO 

150 

STIF0094 

CALL  DC ODE  1 I CODE tltlBCS* 

IBCXt IBCY.1C,IX,IY»NBC) 

STIF0095 

I F 1 I8CS.NE. 1)  bO  TO  110 

STIF0096 

ALF  =  BC 1 IC  t  3 ) 

STIF0097 

BC IX, 1)=BC I X 1 1  I+ALF*  C  ALF*  C  BC I Y , 1  )♦! . ) *2 . *B C I X, 2 )) 

STIF0098 
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B( IX, 2)=-ALF 

STIF0099 

BUY,  1 )  =1. 

STIF0100. 

F (  I X ) = ALF* F (  I Y 1 ♦Ft  IX ) 

STIF0101 

F( IY)=C 

STIF0102, 

KL- IX-MB AND+2 

STIF0103 

KU  = I X+MB AND- 1 

STIF0104 

1FCKL.LT. 1)  KL= 1 

STIF0105 

IF(KU.GT.N2)KU=N2 

STIF0106 

DO  100  K=KL  t  KU 

STIF0107 

IFtK.EC. IX.OR.K.EO. IY)  GO  TO  100 

STIF0108 

IF(K.GT.IY)  GO  TO  90 

STIF0109 

L= I X-K+ 1 

STIF0110 

B(KfL)=B(KfL)+ALF*B(KtL+l) 

STIF0111 

B ( Kt  L+ 1 ) =0 • 

STIF0112 

GO  TO  100 

STIF0113 

90 

L=K- I X+ 1 

STIF0114 

B  (  IX  t  L  )  =B  (  IX,  L)+ALF*B(  IY.L-l) 

STIF0115 

B( IY,L-1)=C. 

STIF0116 

100 

CONTINUE 

STIF0117 

110 

DO  140  N=1 «  2 

STIF0118 

IF(  IFIX(N) .NE.  1)  GO  TO  140 

STIF0119 

IR=  I X+N-l 

STIF0120 

ML= I R-MB AND* 1 

STIF0121 

MU=IR*MBAND-1 

STIF0122 

IFCKL.LT. 11  ML  =  1 

STIF0123 

IF(MU.GT.N2)MU=N2 

STIF0124 

00  130  M  =  ML « KU 

STIF0125 

L=IK-M*1 

STIF0126 

IFCL.LE.l)  GO  TO  120 

STIF0127 

F(M)=F(M)-B(M, L)*BC( IC.N) 

STIF0128 

B ( M, L ) =0 • 

STIF0129 

GO  TO  130 

STIF0130 

120 

L=M- I R+  1 

STIF0131 

F(M)=F(M)-BC IR,L)*BC< ICt N) 

STIF0132 

B ( IR,L)=0. 

STIF0133 

130 

CONTINUE 

STIF0134 

B ( I K  *  1 ) =1 • 

STIF0135 

Ft IR )=BC ( IC, N) 

STIF0136 

140 

CONTINUE 

STIF0137 

150 

CONTINUE 

STIF0138 

RETURN 

STIF0139 

160 

IrfR  l  T  t  (  6 , 1 70  )  M 

STIF0140 

170 

FORKATt 1H010X32HINVALID  ELEMENT  CODE  ELEMENT  NO. 14) 

STIF0141 

STOP 

STIF0142 

END 

STIF0143 
ST  IPO 144 

IBFTC  SYMSL 

SYMSOOOO 

SUBROUT I NE  SYMSOL  ( KKK ) 

SYMS0001 

COMMON/ ADD/  EE( 1C) ,EE1( 10) ,EE2( 10) ,PRR( 10) 

SYMS0002 

COMMON  E,CC,G1,E2, PR.EPR, X21 , Y2 1 , X3 1 , Y3 1, X32 , Y32 , XERR , 

SYMS0003 

1  N2  >  NIELcM»KELt  I  LAW,  MAT  ,  NBC, 

SYMS0004 

2  £:  (  45  0, 22  )  ,  PC ( 30, 3 )  ,  TAB  (  1C1 «  20)  ,  IFIXI2)  , 

SYMS0005 

3  X(45U) ,XCOKD( 225) ,Y( 225) ,IC0DE(225) , 

SYMS0006 

4  FP(450) ,F(450) , 

SYMS0007 

5  11(400) , 12(400) , 13(400) ,NTYPE(400) ,Z(4G0) , 14(400) , 

SYMS0008 

6  SEF( 400) ,SET (400) ,EEP(400) ,EXPL(400) , E YP { 400 1 , EXYP { 400) , 

SYMS0009 

7  NNCOE  *  MBAND 

SYMS0010 

SYMS0011 

NN  =  N2 

SYMS0012 

MM=MBAND 

SYMSOO 1 3 

GO  TO  ( 10,60) ,KKK 

SYMS0014 
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C 

C 

c 


c 

c 

c 


c 

c 

c 


REDUCE  MATRIX 

10  DO  50  N  =  1 *  NN 
DO  40  L  =  2  ,  MM 
C=B(N,L)/B(N,1) 

I =N+L-1 

IF(NN-I)  40,20,20 
20  J=0 

DO  30  K=L , MM 
J  =  J  +  1 

30  B( I ,  J )  =  8  ( I ,J)-C*B(N,K) 

40  B I N , L ) =C 
50  CONTINUE 
GO  TO  130 

REDUCE  VECTOR 

60  DO  80  N= 1 ,NN 
DO  70  L  =  2  * MM 
I =N  +  L- 1 

IFINN-I  )  80,70,70 
70  X  ( I ) -X  (I  ) — B ( N , L ) *X  ( N > 
80  X  (N)=X  (N)/B(N,1) 

BACK  SUBSTITUTION 

N  =  NN 

90  N=N-1 

I  F ( N )  100,  130,100 
100  DO  120  K=2 , MM 
L-N+K- 1 
IF(NN-L) 

110  X  ( N )  =X 
120  CONTINUE 
GO  TO  90 


120,110,110 

(N)-b(N,K)*X  (L) 


130  RETURN 
END 

1BFTC  ST RN 

SUBROUTINE  STRAIN 

COMMON/ ADD/  E E ( 10 ) , EE  1 ( 10 ) , EE2 (  10 ) , PRR ( 10 ) 

COMMON  E,CC,G1,E2,PR,EPR,X21,Y21,X31,Y31,X32,Y32,XERR, 

N2 , NELEM, KEL , I  LAW, MAT , NBC , 

8(450,22) ,BC( 30,3) ,TAB( 101,20) , IFIX(2) , 

X ( 450 ) , XCQRD( 225) ,Y(225) , I  CODE ( 225 ) , 

FP(450),F(450), 

11(400) , 12(400) , I  3(400 ),NTYPE( 400) ,Z(400) ,14(400) , 
SEF(400) , SET (400) ,EEP(400) , E XPL ( 400 ) , E YP ( 400 ) , E X YP ( 400 ) , 
NNODE , MBAND 
130  K=l, NELEM 


1 

2 

3 

4 

5 

6 
7 


DO 

J 1  =  2* 1 1 ( K ) - 
J  2  =  2* 1 1 ( K ) 

J3  =  2* I  2 ( K ) - 1 
J4=2* I  2 ( K ) 
J5=2*I3(K)-1 
J6=2*I3(K) 
CALL  ELEM(K) 
NTY=NTYPE(K) 
E=EE ( NT Y ) 

E 1=EE 1 ( NTY ) 


SYMS0015 
SYMS0016 
SYMS0017 
SYMS0018 
SYMS0019 
SYMS0020 
SYMS0021 
SYMS0022 
SYMS0023 
SYMS0024 
S YMS0025 
S YMS0026 
SYMS0027 
S YMS0028 
S YMS0029 
SYMS0030 
SYMS0031 
SYMS0032 
SYMS0033 
SYMS0034 
SYMS0035 
SYMS0C36 
SYMS0037 
SYMS0038 
SYMS0039 
S YMS0040 
SYMS0041 
SYMS0042 
SYMS0043 
SYMS0044 
SYMS0045 
S YMS0046 
SYMS0047 
SYMS004B 
SYMS0049 
S YMS0050 
SYMS005 1 
SYMS0052 
SYMS0053 
SYMS0054 
S  TRNOOOO 
STRN0001 
STRN0002 
STRN0003 
STRN0004 
STRN0005 
STRN0006 
STRN0007 
STRNOOOB 
S  TRN0009 
STRN0010 
STRN0011 
STKNOO 1 2 
STRNOG 13 
STRN0014 
STRN0015 
STRN0016 
STRN001 7 
STRN001B 
STRN0019 
STRN0020 
STRN0021 
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E2=EE2 (  NT Y ) 

STRN0022 

CC=E1/E 

STRN0023 

PR=PRR(  NTY) 

STRN0024 

IF ( ILAW.EQ.l)  Gl= ( 7«*E/3* )♦*( 1 «. /E2 ) *E 1** < 1 .-l . /E2 ) 

STRN0025 

IF< ILAW.GT.l)  G1=E 1 

STRN0026 

EPR=E/< l.C-PR*PR) 

STRN0027 

I F ( 1 3 ( K ) . E 0 . 0 )  GO  TO  60 

STRN0028 

STRN0029 

****  TRIANGULAR  ELEMENT  CALCULATIONS  **** 

STRN0030 

STRN0031 

10 

A123=X21*Y31-X31*Y21 

STRN0032 

SN=A123/ABSI A123) 

STRN0033 

EXT=<-Y32*X< J1)+Y31*X< J3)-Y21*XI J5) I/A123 

STRN0034 

EYT=I  X32*XI J2)-X31*XI J4)+X21*XI J6) J/AI23 

STRN0035 

EXYT=CX32*XC J1 )-Y32*XC J2)-X31*XI J3)+Y31*XI J4) 

STRN0036 

1 

♦X21*X< J5)-Y21*XI J6) I/AI23 

STRN0037 

EXE  =  EXT-EXPL  <K) 

STRN0038 

EYE=EYT-EYP!K) 

STRN0039 

EXYE=EXYT-EXYPIK) 

STRN0040 

SX=EPR*CfcXE+PR*EYE) 

STRN0041 

SY=EPR* ( EYE+PR*EXE ) 

STRN0042 

SXY=E/I 1.0*PR)*EXYE/2.0 

STRN0043 

SE=SURT<SX**2-SX*SY+SY**2+3.0*SXY**2) 

STRN0044 

CRIT=ABSISE)-SEFCK) 

STRN0045 

IFICRIT ) 40  »  40*  20 

STRN0046 

20 

EET=Sfc/E+EEPlK) 

STRN0047 

CALL  STRSTNIEET,EEPK,SETK*NTY) 

STRN0048 

DEEP=EEPK-EEPIK) 

STRN0049 

30 

EEPIK)=EEPK 

STRN0050 

SET  I K )  =  SfcTK 

STRN0051 

EXPLIK) =DEEP/SE*I SX-SY/2.0)+EXPLIK) 

STRN0052 

LYPIK)=DEEP/SE*ISY-SX/2.G)+EYPIK) 

STRN0053 

EXYP I K ) =3 .0*DEEP/SE*SXY+EXYPI K) 

STRN0054 

ERR=E*DEEP/SE 

STRN0055 

ERR  =  ABS  C  ERR ) 

STRN0056 

GO  TO  50 

STRN0057 

40 

ERR=G.O 

STRN0058 

50 

CONTINUE 

STRN0059 

Q1=E*ZIK)/I 1.0-PR**2)/2.0*SN 

STRN0060 

Q2=E*ZIK)/ll.0+PR)/4.0*SN 

STRN0061 

EXPT=EXPLIK) 

STRN0062 

EYPT=EYP(K) 

STRN0063 

EXYPT=EXYP IK) 

STRN0064 

FPI J1)=-C1*Y32*EXPT-Q1*Y32*PR*EYPT+Q2*X32*EXYPT+FPI Jl) 

STRN0065 

FPC J2)=  Q1*X32*PR*EXPT+U1*X32*EYPT-G2*Y32*EXYPT+FPI J2) 

STRN0066 

FPI J3)=  G1*Y31*EXPT+U1*Y31*PR*EYPT-Q2*X31*EXYPT+FPI J3) 

STRN0067 

FPC  J4)=-Ql*X31*PR*EXPT-01*X31*EYPT+u2*Y31*EXYPT*FPI J4) 

STRN0068 

FPI J5)=-Q1*Y21*EXPT-Q1*Y21*PR*EYPT+U2*X21*EXYPT*FPI J5) 

STRN0069 

FPI J6)=  Ql*X21*PR*EXPT+Gl*X21*EYPT-U2*Y21*EXYPT*FP(J6) 

STRN0070 

GO  TC  110 

STRN0071 

STRN0072 

****  BAR  CALCULATIONS  **** 

STRN0073 

STRN0074 

60 

EET= 1X21*1 XI J3)-XI JI) )+Y21*IXI J4)-XI J2) ) )/Y32**2 

STRN0075 

STRN=ARS I EET-EXPL I K ) ) 

STRN0076 

S IGN= I EET— EXPL IK)) /STRN 

STRN0077 

SE=E*STRN 

STRN0078 

CRI T=SE-SEF I K ) 

STRN0079 

EET=STRN+EEP IK ) 

STRN0080 

IFICRIT)90,9G,70 

STRN0081 

70 

CONTINUE 

STRN0082 

CALL  STRSTNIEET,EEPK,SETK,NTY) 

STRN0083 
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OEEP=EEPK-E£P(K) 

80  EEP ( K ) =EEPK 
SET ( K ) =S€TK 

EXPL(K)=EXPL(K) +S IGN*DEEP 
EH R=E*DEEP/SE 
ERR=ABS ( ERR ) 

GO  TO  100 
90  ERR=0.0 
100  CONTINUE 

EXP  T=EXPL ( K ) 

Q1=E*Z ( K ) / Y  32 
FP(J1)=FP(J1)  “0 1*X  2 1*E  XPT 
FP(J2)=FP(J2 )— Ql*Y21*EXPT 
FP( J3)=FP( J3  )+Ul*X21*EXPT 
FP( J4 )=FP( J4)+gl*Y21»EXPT 
110  IF(ERR-XERR) 130,130,120 
120  XERR=ERR 
KEL  =  K 

130  CONTINUE 
RETURN 
ENO 

$  I BFTC  STRST 

SUBROUTINE  STRSTNl EET , EEPK, SE  TK,NTY ) 

COMMON/ ADO/  EE ( 1 0 ) t EE 1 ( 10 ) t EE2 ( 10) »PRR( 10) 

COMMON  E,CC,Gl,E2,PR,EPk,X21,Y21,X31 ,  Y3 1 , X32 , Y32 , XCRR , 

1  N2,NELEM,KEL, I  LAW, MAT, NBC, 

2  8(450,22),8C(30, 3) ,TA8(  101 ,20  ,  IF  I  X  (2  )  , 

3  X ( 450  ) »  XCORD  (225)*Y(225),  I  CODE  ( 225  )  * 

4  FP ( 4  50 ) ?  F ( 450 ) » 

5  11(400 ,12(400), 13 (400 ),NTYPE( 400) ,Z( 400) « 14 (400 9 » 

6  SEFI400) ,SEF(400),EEP( 4C0) , EXPL ( 400 ) , EYP( 400 ) , EXYP ( 400 ) , 

7  NNOOE,MBANO 

GO  TU ( 10,50, 60 ) , IL AW 
10  J=5.0*EET/CC+1.0 
NT  2~ 2*NTY 
NT1=NT2-1 
I  F  (  J  -  1 0 1  )  2  0 , 30 , 30 

20  EEPK*TAB ( J , NT2 ) ♦  ( T AB ( J  + 1 , NT2 ) ~T AB < J , NT  2 ) >* ( EE T-T AB ( J ,NT  1 ) ) / 
1(TAB(J+1,NT1)-TA8(J,NT1) ) 

GO  TO  40 

30  CALL  PNEW1 ( tEPK , EET ,E,G1,€2) 

40  SETK=G1*EEPK**< 1.0/E2) 

RETURN 

50  5ETK=E*EET/< l.+(ABS(E*EET/Gl) )**E2)**(1./E2) 

EEPK  =  tET -SE  T  K / E 
RETURN 
60  EC*G1/E 

IF(FfcT-EC)70,70,S0 
70  EEPK=0« 

SETK=E*EET 

RETURN 

80  EEPK=EET-EC 

SETK=G1+E2*EEPK 
m  RETURN 
ENO 

$  I  BFTC  OTPUT 

C 12345  OUTPUT  SUBROUTINE 

SUBROUTINE  CUT P T 

COMMON/ ADO/  E£<  10)  ,EE1(  10)  ,EE2<  10)  ,PRRUO) 

COMMON  E,CC,G1 , E2,PR,EPR, X21 , Y21,X31 , Y3 1 , X32 , Y32 , XERR , 

1  NN2 ,NELEM , KEL, I  LA W  j  M AT ,NBC, 

2  B(450,22),BC(30,3) ,TABt 101,20) , IF  IX (2), 


STRN0084 
STRN0085 
STRN0086 
STRN0087 
STRN0088 
STRN0089 
STRN0090 
STRN0091 
STRN0092 
STRN0093 
STRN0094 
STRN0095 
STRN0096 
STRN0097 
STKN0098 
STRN0099 
STRN0100 
STRN010 1 
STRN0102 
STRN0103 
STRN0104 
STKSOOOO 
STRS0001 
STRS0002 
STRS0003 
ST  KS0004 
STRS0005 
STRS0006 
STKS0007 
STRS0008 
STRS0009 
STKS0010 
STRSOO 1 1 
STKS0012 
STRS0013 
STKS0014 
STRSOO 1 5 
STKSOO  L  6 
STRS0017 
STRS0018 
STRS0019 
S  TRS0020 
STRS0021 
STRS0022 
STRS0023 
S  TKS002  4 
STRS0025 
STRS0026 
STKS0027 
STRS002B 
S  TKS0029 
STKS0030 
STRS0031 
STRS0032 
STRS0033 
OTPUOOOO 
0TPU0001 
OTPU0002 
0TPU0003 
QTPU0004 
0TPU00G5 
0TPU0006 
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3  X ( 450 ) 9  XCORD (225),Y(225) ,ICODE(225) , 

4  FP(450) ,F(450) , 

5  11(400) 9 12(400) , 13(400) » NT YPE( 400) ,Z(400 

6  SEF(400) 9SET(400),EEP(400) , EXPL ( 400 ) , E YP 

7  NNODE  9  MBAND 
L  =  0 

DO  70  K  =  1 *  NELEM 

♦  ♦♦♦  TRIANGULAR  ELEMENT  CALCULATIONS  ♦  ♦♦♦ 

NTY  =  NTY  PE ( K ) 

E=EE ( NT Y ) 

PR=  PRR(NTY) 

EPR=E/( l.-PR*PR) 

I F ( I3(K).EQ.O)  GO  TO  70 
10  CALL  EL  EM  K  ) 

A123=X21*Y31-X31*Y21 

J1=2*I1(K)-1 

J2=2*I1(K) 

J3  =  2* I 2 ( K ) — 1 
J4=2*I2(K) 

J5  =  2*  1 3 ( K )  —  1 
J6  =  2* I  3 ( K ) 

EXT=(-Y32*X( J1 )+Y31*X( J3)-Y21*X( J5) )/A123 
EYT= (  X32*X( J2)-X31*X( J4)  +  X2UX( J6 ) )/A123 
EXYT=(X32*X( J1 )-Y32*X( J2)-X31*X( J3)+Y31*X( J4) 

I  +X21*X( J5)-Y21*X( J6) J/A123 
EXE=tXT-EXPL(K) 

EYE=EYT-E YP ( K ) 

EXYb^EXYT-EXYP(K) 

SX=EPR* ( EXE+PR*EYE) 

SY  =  EPR* ( EY  E+PR*EXE ) 

SXY=l/( 1.0+PR)*EXYE/2.0 

PE2  =  SQRT  (  (  .  5*  (  SX-SY  )  )*^2  +  SXY*»2) 

PHI=.5^ATAN2( (-2.0*SXY) f (SX-SY) ) ♦57.29578 
PH2  =  • 5* AT  AN2 ( (-EXYT) f (EXT-EYT)  )♦ 57. 2957  8 
PS  1  = . 5* ( SX  +  SY ) 

SIGE1=PS1+Pt 2 
S I GE2=PS 1-PE2 
PST1  =  .5M  EXT  +  EYT ) 

PE T2  =  SQRT (  (  •  5*  (  EXT-EY  T  )  )  ♦♦2  +  EX YT^2/4. 0 ) 

STRE1=PST1+P£T2 

STRE2=PST1-PET2 

PET 2  =  2 • O^PET  2 

N 1  =  1 1  (  K  ) 

N2=I2(K) 

N3=  I  3 ( K ) 

XC=XCnRD(Nl) ♦  ( X21+X31 )/3.0 
Y  C  =  Y (Ml )  ♦  ( Y  2  1  ♦  Y  3 1 ) /3.0 

EPE=2 .♦ SORT ( (EXPL(K)+*2  +  EXPL(K)  ♦EYP(K)+EYP(K)«2 
1) 

SE=SuRT ( SX*«2-SX*SY  +  SY^2  +  3.*SXY **Z) 

L  =  L  ♦  1 

I  F ( MUD ( L* 1 9  14) )  30,20  9  30 
20  WRITE  ( 6  9  40 ) 

30  WR I T  c ( 6  9 50 )KfXC9YC9SXf  SY,SXYf  S I GE1 » S I GE2 , P E2 » PH  I 
lYTfSTRElf STKE29PET2 

40  FORMAT (9H1EL.  NO. / 5X1 1MC00RDINATES2GX33HS  T  R  E 

II  N  S  /8H  PHI 7X1HX8X1HY6X9H  TAU-XX6X9H 

2Y8X7HMAX IMUM8X7HM I N I MUM6X9HMAX  SHEAR  ) 

50  FORMAT* lH0I7f0PF8.39F9.39lP6E15.4/lH  0PF7.2,F8.2 
WK I T b ( 6 9  60 )  SE , EPt 


OTPU0007 

OTPU0008 

)  ,  14(400)  9  OTPU0009 

(400) , EX YP( 400) ,  OTPUOOIO 
OTPUOOl 1 
OTPUOO 12 
0TPU0013 
0TPU0014 
OTPUOO 15 
OTPU0016 
OTPU0017 
OTPUOO 1 8 
OTPUOO 19 
OTPU0020 
OTPU0021 
OTPU0022 
OTPU0023 
OTPU0024 
0TPU0025 
OTPU0026 
OTPU0027 
OTPU0028 
OTPU0029 
OTPU0030 
OTPU0031 
0TPU0032 
OTPU0033 
OTPU0034 
OTPU0035 
OTPU0036 
OTPU0037 
OTPU0038 
OTPU0039 
OTPU0040 
OTPU0041 
OTPU0042 
OTPU0043 
OTPU0044 
OTPU0045 
OTPU0046 
OTPU0047 
OTPU0048 
OTPU0049 
OTPU0050 
OTPU0051 
OTPU0052 
OTPU0053 
0TPU0054 
OTPU0055 

+EXYP ( K )*+2/4. ) /3.0TPU0056 
OTPU0057 
OTPU0058 
OTPU0059 
OTPU0060 
OTPU0061 

, PH2  9EXT9EYT9  EXOTPU0062 
OTPU0063 

SSES/STRA  OTPU0064 

TAU-YY6X9H  T  AU-XOT  PU0065 

OTPUOO66 

, 9X1 P6E 15 • 4 )  OTPU0067 

OTPUOO68. 
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60  FORMAT ( 1H 
1N=E12.5) 
70  CONTINUE 


23H*****  EFFECTIVE  STRESS=E 1 2 . 5, 23H*****  EFFECTIVE 


C 

C 

C 


****  BAR  CALCULATIONS  **** 

J=0 

DO  120  K=1 , NELEM 
NTY=NTYPE ( K ) 

E  =  EE (  NTY ) 

PR=PRR (  NTY) 

J 1  =  2* 1 1 ( K )  —  1 
J  2  =  2* 1 1 ( K ) 

J  3  =  2* 1 2 ( K ) - 1 
J4=2*I2(K) 

I  F  ( I  3 ( K  )  .NE • 0 )  GO  TO  120 
80  CALL  ELEM(K) 

EET=(X21*(X( J3)-X( Jl) )+Y21*(X(J4)-X(J2)  )  )/Y32**2 
SE=E*(EET— EXPL(K) ) 

SEMZK=SE*Z(K) 

K1=I1(K) 

K2=I 2  ( K  ) 

I  F  ( J ) 100t  90f 100 
90  WRITE  (6,130) 

J  =  1 

100  CONTINUE 

110  WRITE(6,140)K,Kl,K2fSE,EtT, SEMZK 
120  CONTINUE 

130  FORMAT ( 9H0B A R  NO.  6X 10HN0DE  NOS. 

1 13HMEMBER  FORCES) 

140  FORMAT ( 1H0  3  1 8 , 2E 1 5 . 5 , 2X , E 15 . 5 ) 

RETURN 
END 

$DAT  A 

NONLINEAR  TRUSS  PROC.  ASCE  DEC.  1965 


8X7HSTRESS  8X7HSTR A  I N , 4X , 


STRAI 0TPU0069 
0TPU0070 
OTPU0071 
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0TPU0075 
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0TPU0089 
0TPU0090 
0TPU0091 
0TPU0092 
0TPU0093 
0TPU0094 
0TPU0095 
0TPU0096 
OTPU0097 
OTPU0098 
0TPU0099 
OTPUOIOO 
0TPU0101 


ONE  MATERIAL  RAMBERG  OSGOOD  LAW 


6  10  1 
0.10000000E+05 


10  20 

1  1 

2 

3 

4  1 

5 

6 

1  4 

2  4 

3  1 

4  1 

5  2 

6  5 

7  5 

8  2 

9  2 

10  3 

3 

NONLINEAR 
6  10  1 
0 • 10000000E+05 
0 • 10000000E+05 
10  20  1 

1  1  1 

2 


I  1 

40.50000 
1  10 

30. 

60. 


22 

0.3 

1 


30 


.01 


7.0 


5 

2 

5 

2 

5 

6 
3 
6 
3 
6 


0 
0 
0 
0 
0 
0 
0 
0 
0 
0 

-10. 

TRUSS  PROC 
G  2 
65.00000 
45.00000 
1  10 


30. 

60. 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


40. 

40. 

40. 


.25 

.20 

.20 

.25 

.20 

.25 

.20 

.20 

.25 

.20 


ASCE  DEC.  1965 
22  30 

C.3 


0. 

1 


TWO  MATERIALS  RAMBERG  OSGOOD  LAW 

8.0 

8.0 


.01 


30. 


45 
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3 

60. 

4 

1 

1 

40. 

5 

30. 

40. 

6 

60. 

40. 

1 

4 

5 

0  1 

.25 

2 

4 

2 

0  2 

.20 

3 

1 

5 

0  2 

.20 

4 

1 

2 

0  l 

.25 

5 

2 

5 

0  1 

.20 

6 

5 

6 

0  1 

.25 

7 

5 

3 

0  2 

.20 

8 

2 

6 

0  2 

.20 

9 

2 

3 

0  1 

.25 

10 

3 

6 

0  1 

.20 

3 

-10. 

MIT 

SHEAR 

LAG 

PROBLEM 

CONFIGURATION 

49 

78 

2 

1 

22  30 

10200000E+08 

52000. 

0.3  5 

10 

20 

1 

1 

1  0.01 

1 

1 

1 

G.G 

2 

1 

0.5 

3 

1 

1.0 

4 

1 

1.5 

5 

1 

2.0 

6 

1 

0.0 

0.5 

7 

0.5 

0.5 

8 

1.0 

0.5 

9 

1.5 

0.5 

10 

2.0 

0.5 

11 

1 

3.0 

12 

1 

4.0 

13 

1 

5.0 

14 

I 

0.0 

1.0 

15 

0.5 

1.0 

16 

1 .0 

1.0 

17 

1.5 

1.0 

18 

2.0 

1.0 

19 

3.0 

1.0 

20 

4.0 

1.0 

21 

5.0 

1.0 

22 

1 

0.0 

2.0 

23 

0.5 

2.0 

24 

1.0 

2.0 

25 

2.C 

2.0 

26 

3.0 

2.0 

27 

4.0 

2.0 

28 

5.0 

2.0 

29 

1 

0.0 

3.0 

30 

0.5 
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31 

1  .0 
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32 
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3.0 

33 
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3.0 

34 
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35 
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3.0 

36 

1 
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4.0 

37 

0  •  t) 

4.0 

38 

1.0 

4.0 

39 

2.0 

4.0 

40 

3.0 

4.0 

41 

4.0 

4.0 

42 

5.0 

4.0 

43 

1 

0.0 

5.0 

46 
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5.0 
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1 

2 
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2 

2 

6 

7 

1 
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3 

2 

3 

7 

1 
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4 

3 

7 

8 

1 

•  08 

5 

3 

4 

8 

1 

.08 

6 

4 

8 

9 

1 

.08 

7 

4 

5 

9 

1 

•  08 

8 

5 

9 

10 

1 

.08 

9 

5 

10 

11 

1 

.08 

10 

10 

11 

19 

1 

.08 

11 

11 

12 

19 

1 

•  08 

12 

12 

19 

20 

1 

.06 

13 

12 

13 

20 

1 

.08 

14 

13 

20 

21 

1 

.08 

15 

6 

7 

14 

1 

.08 

16 

7 

14 

15 

1 

.08 

17 

7 

8 

15 

1 

•  08 

18 

8 

15 

16 

1 

.  08 

19 

8 

9 

16 

1 

.08 

20 

9 

16 

17 

1 

.08 

21 

9 

10 
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1 

•  08 

22 

10 
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18 

1 

.08 
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10 
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23 
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19 
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1 
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20 

26 

27 

1 
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35 

20 

21 

27 

1 

•  08 

36 

21 

27 

28 

1 

.08 

37 

22 

23 

29 

1 

.08 

38 

23 

29 

30 

1 
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39 

23 

24 

30 

1 
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40 

24 

30 

31 

1 

.08 

41 

24 

25 

31 

1 

.08 

42 

25 

31 

32 

1 
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43 

25 

26 

32 

1 

.06 

44 

26 

32 
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1 

•  08 

45 

26 

27 

33 

1 

.08 

46 

2  7 

33 

34 

1 
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27 

28 

34 

1 

.06 

48 

28 

34 

35 

i 

.08 

49 

29 

30 
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